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We present a fully dynamical model to study nonequilibrium effects in both the chiral 
and the deconfinement phase transition. The sigma field and the Polyakov loop as the 
corresponding order parameters are propagated by Langevin equations of motion. The locally 
thermalized background is provided by a fluid of quarks and antiquarks. Allowing for an 
exchange of energy and momentum through dissipative and stochastic processes we ensure 
that the total energy of the coupled system remains conserved. We study its relaxational 
dynamics in different quench scenarios and are able to observe critical slowing down as well as 
the enhancement of long wavelength modes at the critical point. During the fluid dynamical 
expansion of a hot plasma fireball typical nonequilibrium effects like supercooling and domain 
formation occur when the system evolves through the first order phase transition. 

PACS numbers: 25.75.-q, 47.75.+f, 11.30.Qc, 24.60.Ky, 25.75.Nq 
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I. INTRODUCTION 

One of the major goals of heavy-ion physics is to gain firm knowledge about the different phases 
of strongly interacting matter and the transitions between them. At high temperatures, the chiral 
symmetry that is broken in the QCD vacuum, is expected to be restored. Furthermore, quarks 
and gluons might become the relevant degrees of freedom after the transition from hadrons to a 
deconfined state of matter. From lattice QCD calculations we know that at vanishing baryochem- 
ical potential the chiral transition is an analytic crossover [HE]- As recent studies show, there 
seems to be no direct relation between the chiral restoration and the onset of deconfinement j3j- 
[5]. A multitude of effective models have been used to investigate the QCD phase diagram [6HH] 
indicating a first order phase transition at large ending at a critical point (CP). In equilib- 
rium this CP exhibits a divergence in the correlation length of the order parameter. The problem 
how to locate the CP in the T-/X£-plane in heavy-ion collisions was addressed by Stephanov, Ra- 
jagopal and Shuryak \15\ [TC] for equilibrated systems. They proposed to search for divergences in 
event-by-event fluctuations of quantities like transverse momentum or particle multiplicity, since 
the strength of these fluctuations is directly related to the correlation length of the sigma field, 
the order parameter for the chiral phase transition. However, a priori one cannot expect to reach 
thermodynamic equilibrium during such a heavy-ion collision, especially near the phase transition 
where relaxation times become large in comparison to the rapid dynamics of the expanding fireball. 
Therefore, several additional effects have to be taken into account. Not only is the growth of the 
correlation length naturally limited by the finite size of the system, but also critical slowing down 
of the dynamics in the vicinity of the critical point will weaken the expected signals [T7] . One may 
try to overcome this difficulty by looking for quantities that are accessible in experiment and more 
sensitive to the correlation length. Higher cumulants and combinations of them like the kurtosis of 
conserved quantities such as net baryon number or net charge fulfill this requirement [T9"] . The 
fast dynamical evolution of matter through a first order phase transition may exhibit interesting 
phenomena like supercooling followed by a decay through spinodal decomposition |20H25| or nu- 
cleation [26J . Also an enhancement of soft pions from the decay of a disoriented chiral condensate 
(DCC) was proposed ELS 3> SI gnal for a nonequilibrium chiral phase transition [27-29J. 

Only via thorough theoretical understanding of the processes and effects in dynamical systems 
undergoing phase transitions one may provide realistic predictions for upcoming experiments with 
relativistic heavy-ion beams at RHIC [30J, CERN SPS [31], FAIR [32] and NICA [33]. 

In [8] , the thermodynamic properties of a Polyakov loop extended Nambu- Jona-Lasinio (PNJL) 
model were studied and compared with lattice QCD results. The addition of the Polyakov loop to 
these models improved the behavior of thermodynamic bulk quantities in comparison to the lattice 
QCD data. A similar extension to the sigma model with constituent quarks was presented in [9] 
where a coincidence of the deconfinement and chiral symmetry transition even at finite \ib was 
found. An investigation of the phase structure of that model beyond mean-field has been performed 
in [10] including quantum fluctuations within the functional renormalization group method. This 
improvement has lead to a quantitative shift in the position of the CP. Furthermore, the authors 
calculated net-quark number density fluctuations as well as ratios of cumulants as an important 
means to identify the location of the deconfinement and chiral phase transitions. The relevance 
of the fermion vacuum loop for such models has been investigated in [11] where it was shown 
that this term has crucial influence on the phase structure and on physical observables such as 
net quark number fluctuations. In [12] an extension of this model to (2+1) flavors together with 
a comparison with lattice QCD data have been presented. An alternative model using a dilaton 
field representing a scalar glueball condensate was subject to studies in [13] . An important step 
towards the understanding of nonequilibrium effects at the chiral phase transition was done in |14] , 
where it was shown that the inclusion of spinodal instabilities yields an enhancement of density 
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fluctuations along the first order transition line. This is in contrast to the usual understanding of 
thermalized systems where such fluctuations are expected to diverge or grow large only at the CP. 

For a better understanding of the processes during a heavy-ion collision nonequilibrium effects 
have to be taken into account in the framework of a fully dynamical model. A promising ansatz 
to study the QCD phase transition in such a way is given by a chiral fluid dynamics model that 
has been developed and extended over more than ten years now [34H38] . In |34| . the authors 
introduced a model coupling a relativistic ideal fluid of quarks to the linear sigma model and a 
scalar glueball condensate. Later in |35j this model was used to study the 3 + 1 dimensional fluid 
dynamic expansion of a plasma droplet coupled to the out-of-equilibrium evolution of the long- 
wavelength modes of the chiral condensate. There the initial fluctuations of the sigma field were 
propagated deterministically. A self-consistent derivation of the coupled dynamics of fields and 
fluid together with the local thermodynamic properties of the heat bath was then presented in [36] . 
This step is a crucial improvement as the previous models which propagated the chiral fields using 
a classical equation of motion neglected relaxational and stochastic processes. As as result, the 
fields would never relax to their equilibrium state for a given temperature but continue oscillating. 
In [36], the order parameter of chiral symmetry is propagated by a Langevin equation to include 
damping and noise in the heat bath of quarks. The quark fluid is propagated according to the 
equations of ideal hydrodynamics. The relaxational dynamics of the field and the treatment of 
the finite size of the heat bath were investigated in [37] and the first numerical studies of the full 
expanding system were then presented in |38| . A coupling of the linear sigma model to viscous 
hydrodynamics has been studied in [39J. 

Dynamical models for the Polyakov loop have been proposed in [40H42] . In |40| . the authors 
developed a model for particle production near the deconfinement phase transition due to oscil- 
lations in a Polyakov loop condensate. These oscillations were included on the basis of a kinetic 
term in an effective field theory for a Polyakov loop that is biquadratically coupled to a mesonic 
field. This idea was further used in |41| . where the authors showed how explosive behavior at the 
QCD phase transition might be produced by the decay of such a condensate of Polyakov loops. 
Based on work in |43[ 01] , a Langevin equation for the deconfinement order parameter for pure 
SU(2) gauge theory was developed in [42] . It was shown that the dissipative interaction with the 
medium plays a significant role in the determination of physical time scales. 

In the present work, we connect these ideas of a dynamical Polyakov loop with the model 
of chiral fluid dynamics to take into account effects of the deconfinement phase transition. The 
newly included deconfinement order parameter is treated as an effective field and propagated 
by a phenomenological Langevin equation. The suggested model is able to provide a dynamic 
description of two phase transitions in the background of a fluid dynamically expanding heat bath 
of quarks. This setup resembles the situation after the collision of two heavy nuclei. In |36j the 
proper nonequilibrium dynamics for the sigma field and the quark fluid have been derived self- 
consistently for a model without Polyakov loop. The resulting Langevin equation for the sigma 
field includes two important effects: the damping of the chiral field due to the interaction with 
the heat bath and the back reaction of the heat bath on the sigma field by stochastic noise. For 
the dynamics of the Polyakov loop we follow the same idea but on a phenomenological basis as 
it is currently not possible to derive the dynamics of the Polyakov loop from a consistent field 
theoretical approach like it has been achieved for the sigma field. With the considered Polyakov 
loop extended linear sigma model we are able to describe characteristic phenomena at the phase 
boundary within a dynamical setup. 

This article is structured as follows. In section[ll]we present the model of Polyakov loop extended 
chiral fluid dynamics, where the sigma field and Polyakov loop are coupled to a fluid dynamically 
propagated heat bath of quarks. After that we consider two different numerical implementations. 
Section III presents results of relaxational dynamics in a box for several quench scenarios. In 
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section IV we study the evolution of a freely expanding hot plasma undergoing a phase transition. 



Finally, we present a short summary and outlook in section \V\ 

II. CHIRAL FLUID DYNAMICS WITH A POLYAKOV LOOP 

A. General remarks 

Our model extends existing studies of chiral fluid dynamics [34-38J with the Polyakov loop to 
take into account effects of both the chiral and the deconfinement phase transition. The restoration 
of SU(./Vj)l x SU(iVf )r chiral symmetry in the high temperature phase can be characterized by the 
melting of the (qq) condensate or the sigma field, respectively. For the transition to deconfinement, 
the Polyakov loop £ is the quantity of interest. It is defined as the expectation value of the color 
trace of a thermal Wilson loop: 



t= Wc (tr c V)p, £=—(tr c V j r)p, (1) 
where the operator V is defined as 

V = Pexp(ig B I drAo] . (2) 



\ Jo / 

Here, P stands for path ordering, Aq denotes the temporal component of the Euclidean color gauge 
field, g s the strong coupling constant and ft = 1/T is the inverse temperature. The expectation 
value of £ is related to the free energy F q (x) of an infinitely heavy test quark at spatial position x 
via 

(£(x)} = e-WW . (3) 

In the confined phase, this free energy diverges and therefore (£) vanishes whereas it takes some 
finite value in the deconfined phase. In the limit of infinitely heavy quarks, QCD is invariant under 
Z(N C ) center symmetry transformations of the SU(N C ) color gauge group. The Polyakov loop 
transforms as 

£ -> zl , z G Z(N C ) , (4) 

with the consequence that the confined phase is center symmetric while this symmetry is sponta- 
neously broken in the deconfined phase. In the presence of dynamical quarks, the Polyakov loop 
is always non-zero as the free energy of a test quark does not diverge anymore. Nevertheless, (£) 
may still serve to characterize the transition between the two phases. 



B. The model 

For our studies we use the Polyakov loop extended quark meson model [SJ. The Lagrangian 
reads 

C = q [i (Vfy - ig s l°A ) -g(a + i lh r ■ if)] q + i (<V) 2 + I (d M vr) 2 - U (a, if) - U(£, £) , (5) 

where q = (u,d) is the constituent quark field, so Nf = 2, and a the mesonic sigma field. For 
our simulations around the phase transition we utilize a fixed strong coupling constant of a s = 
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gs/(47r) = 0.3. As we are only interested in the behavior of the order parameters, we neglect 
fluctuations of the pionic degrees of freedom and keep their values fixed at the vanishing expectation 
value 7? = (tt) = for all times. The potential for the sigma field reads 



U(a) = ^{a 2 -u 2 ) 2 - h q a-U . (6) 



The chiral symmetry of the Lagrangian ([5j) is explicitly broken by the term h q in the potential ([6j) 
taking into account the finite quark masses. The parameters in ^ are chosen such that chiral 
symmetry is spontaneously broken in the vacuum, where (a) = f n = 93 MeV, the pion decay 
constant. The explicit symmetry breaking term is h q = f n m 2 with the pion mass m n = 138 MeV. 
These requirements lead to v 2 = f 2 — m 2 /X 2 . Choosing A 2 = 19.7 yields a realistic vacuum sigma 
mass m 2 = 2X 2 f 2 + m 2 « 600 MeV. The constant term Uq = m^./(4A 2 ) — f 2 m 2 is chosen such that 
the potential energy vanishes in the ground state. The quark-meson coupling constant is fixed by 
the requirement to reproduce the constituent quark mass in vacuum, g = m q / = 3.3. 

The temperature dependent Polyakov loop potential is chosen in a polynomial form [8, 9, 45J. 

£ M = -^2 (\H 2 + \i\ 2 ) " | (i 3 + ?) + ^ (\H 2 + |f) 2 ■ (7) 

At non- vanishing ^5, this parametrization yields an unphysical behavior in the Polyakov loop 
susceptibilities which become negative in a broad temperature range [46 i. One can cure this 
problem by augmenting this effective potential with a logarithmic term to account for the Haar 
measure in the group integral [46--48J . As we will restrict our numerical studies to the case of zero 
baryochemical potential, we can ignore this issue for the present investigation. 

The coefficients in ([7]) are fixed to reproduce thermodynamic results from lattice QCD simula- 
tions in the pure gauge sector. We use parametrizations proposed in [H l4"9"H5"2] 

6 2 (T)=a + ai (|) + a 2 (|) 2 + a3(|) 3 (8) 

with ao = 6.75, a\ = —1.95, 02 = 2.625, 03 = —7.44 and two temperature independent coefficients 
63 = 0.75 and 64 = 7.5. The potential ([7]) has a first-order phase transition at a critical temperature 
of To = 270 MeV. However, in the presence of dynamical quarks this transition temperature 
gets lower due to the increased number of degrees of freedom and reduced dynamical scale. The 
proper Nf- and //^-dependence of To has been investigated in [53]. For two flavors and vanishing 
baryochemical potential, To reduces to a value of 208 MeV. 

The effective thermodynamic potential that we need for describing the dynamics of a and I and 
for the local equilibrium properties of the quarks is given by 

V e $ = -ylnZ = U + U + n qq . (9) 

The partition function Z of our system can be written as a path integral over the quarks, antiquarks, 
mesons and the temporal component of the color gauge field. Integrating out the quark degrees 
of freedom, which will constitute the heat bath, we obtain the grand canonical potential £l qq . At 
li B = 0, £ = £ and in the mean-field approximation it is [9] 



4NfT J ^3 In [l + Me-? E + 3£e~ 2 ^ E + e" 3 ^] . (10) 



Here E = y 1 'p 2 + g 2 a 2 , the energy of the quarks, their mass being generated dynamically by the 
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sigma field. The integrand contains contributions from one-, two- and three-quark states, the 
first two being proportional to I. This means that for vanishing value of the Polyakov loop only 
three-quark states contribute, while the amount of one- and two-quarks gets larger with growing I. 
This is called "statistical confinement" [S3]. In (10), we omit the zero-temperature contribution to 
Qqq which can partly be renormalized into the parameters A 2 and u 2 , leaving a logarithmic term 
depending on the renormalization scale and the effective quark mass. This term may have crucial 
influence on the phase structure of the model 155]. However, as the mean-field approximation 
provides us with the desired phase transition already [6], we neglect this term and its effects in 
the following. In order to simplify the calculations and for a first qualitative study we follow the 
same strategy as in [35H37] . Varying the coupling strength g one can tune the characteristic shape 
of the effective potential V e s at [is = and by that the type of transition: For g = 4.7 we see 
two degenerate minima (a = 9 MeV,^ = 0.40) and {a = 81 MeV,£ = 0.22) at the transition 
temperature of T c = 172.9 MeV, see fig. [la] While for g = 3.52 we have only one single minimum 
(a = 49 MeV, I = 0.43) at T c = 180.5 MeV, where the potential is very broad and flat as shown in 
fig. lb This resembles a CP. Note that in principle one has to choose g such that the product ga 
in vacuum reproduces the constituent quark mass, which leads to a value of g = 3.3. 
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FIG. 1. la Effective potential for g — 4.7, corresponding to a first order phase transition at T c 
[Tb] Effective potential for g = 3.52, corresponding to a CP scenario at T c = 180.5 MeV. 



172.9 MeV. 



C. The equations of motion 

In [3B] we have derived the coupled dynamics of the sigma field and the quark fluid self- 
consistently with the two particle irreducible (2PI) effective action for an analogous model without 
Polyakov loop. The description of nonequilibrium processes can be achieved via Langevin equa- 
tions. This has been extensively done in a quantum field theoretical framework for </> 4 theory 
[S6H59] . gauge theories [60] [61] and O(N) chiral models [62]. Here, a splitting between the long- 
and short-wavelength modes of the sigma field was assumed and the relaxational dynamics of the 
soft modes in the heat bath of hard modes was derived within the influence functional formalism. 
Utilizing a chiral model with constituent quarks, we assumed a different splitting. Taking the 
quarks as the environmental degrees of freedom and the sigma as the relevant degrees one can 
calculate the 2PI effective action by integrating over the Keldysh contour. Out of that we were 
able to derive a Langevin equation of motion containing friction and noisd^] 

1 In this form the equation of motion is not Lorentz invariant. The problem could be cured by replacing r\ a dto —¥ 
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d^a + Vf7 (T)d t a + 



dV c 



da 



(ii) 



The damping coefficient r) a is temperature dependent and responsible for the transfer of energy- 
momentum to the heat bath. It arises from the a qq process and is non-vanishing wherever this 



decay is kinematically possible. In |36j we derived its explicit form in the 
sufficient as we are interested mainly in the fluctuations of the soft modes 
given by 



limit which is 



tor m r 



> 2m q it is 



12 5 2 



7T 
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m. 



mi 



(12) 



with m a being the pole mass of the sigma field. While the screening mass, defined as the curvature 
of the effective potential at the equilibrium value, vanishes at the CP, this is not in general true 
for the pole mass. However, we use 



m„ 



d 2 V eS 



da 2 



(13) 



as an approximation to the proper pole mass for equation (12). This ansatz is justified as the pole 



mass has a minimum at the transition temperature for a crossover scenario just like the screening 
mass. A consistent calculation is possible by considering the dispersion relations for different types 
of excitations. 

For the damping below the phase transition we choose rj a = 2.2 /fm |28j to account for the 
a f-> 7T7T reaction whenever it is kinematically allowed. Fig. [2] shows rj a as a function of temperature 
for both transition scenarios. Due to its perturbative derivation, the damping coefficient is rather 
high, especially for the first order transition. The only region where it vanishes is around T c for 
the CP, where the sigma becomes light enough so that m a < 2m q , 2m n . The stochastic noise field 
in the equation of motion (11) has a vanishing expectation value = 0. To complete this part 
we derive a dissipation-fluctuation-relation that connects the damping coefficient to the correlator 
of the stochastic noise field and ensures the relaxation to the proper equilibrium state. 



(Ut, S)Ut', #)) = ^8{t - t')6(x - ^)m aVtJ coth (^) . (14) 



The temperature, the mass of the sigma field, the damping and the noise field in ( 14 ) are local 
quantities which may vary between different positions on the spatial grid. 

The derivation of the equation of motion for the Polyakov loop I is more problematic because 
t is defined in the Euclidean space and it is not completely clear how to propagate it in real time. 
In [3D] a kinetic term for the Polyakov loop has been included in the Lagrangian: 



C + -^\d^\ 2 T 2 . (15) 

It was proposed that production of pions near the phase transition is driven by oscillations of the 
Polyakov loop. Only in the region around T c the Polyakov loop is light, allowing large fluctuations 
and thus particle production. As the Polyakov loop operator is a phase in color space and therefore 
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FIG. 2. Damping coefficient for the sigma field as a function of temperature, for a CP and a first order 
transition scenario. 



£ a pure number, dimensions in an effective Lagrangian can only be made up by powers of the 
temperature T. However, the temperature in our model is space and time dependent, so the 
Euler-Lagrange equation for £ would contain derivative terms of T that might lead to unphysical 
behavior. 

For our present model we follow a different strategy formulated in [63]. If the order parameter 
is out of equilibrium, then relaxational processes towards its thermodynamic equilibrium value will 
occur. The velocity of relaxation is then proportional to the derivative of the effective potential 
with respect to the order parameter. Analogous to the sigma field we include fluctuations of the 
Polyakov loop with a stochastic noise term for which = 0: 



m d t £T 2 + = ^ , (16) 



Mt, *Mt', x'))T 2 = y5(t - t')8{x - x')2 Ve T . (17) 



Equation (17) assumes Gaussian and Markovian noise as it was used in [12]. In contrast to the 
sigma field it is not clear how to derive a damping for the Polyakov loop in a field theoretical 
manner as the fermionic part of the Lagrangian contains interaction with the Aq color gauge field 
but not directly with £. Therefore we chose a 'reasonable' value of rji = 5/ fm, although the results 
are not sensitive to it. We will further comment on this later. 



D. Propagation of the quark fluid 



The propagation of the quark fluid is governed by the equations of ideal relativistic fluid dy- 
namics 

d M (Tf + 7T + If) = . (18) 

We require energy-momentum conservation for the system as a whole, its energy- momentum tensor 
being split into quark, sigma and Polyakov loop parts. In [36] we were able to derive self-consistently 
an approximate expression for the quark and sigma contribution. For the quarks we obtained the 
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energy-momentum tensor of an ideal fluid, the contribution of the sigma field can be expressed as 

d^=(-°^-T h d t< r\ff'*. (19) 
Assuming the same structure and effects for the Polyakov loop contribution, we end up with 

W = (-^jjr " ^T 2 ) ■ (2°) 

As explained in [36 , these terms cannot account for the average energy transfer from the heat 
bath to the fields caused by the stochastic noise fields £ CT and ^. The equation of state p = p(e) is 
obtained and tabulated from the thermodynamic relations 

p(a,£,T) = -n g g(a,e,T) , (21) 

e(a,i,T)=T dpia d ^ T) -p(a,£,T). (22) 

A simple relation for p{e) cannot be obtained as a and i are propagated explicitly and can be out 
of equilibrium and the pressure as well as the energy density of the quarks depend explicitly on 
the local values of these fields. 



III. EQUILIBRATION IN A BOX 

In this section we study the relaxational behavior of the order parameters in a box of finite size 
after several temperature quenches. The aim is to give estimates and compare relaxation times 
near and far away from the transition point for a first order and a CP scenario. Furthermore we 
investigate fluctuations during and after the equilibration process. 



A. Numerical implementation 

We solve the equations of motion for the fields and the quark fluid on a fixed spatial cube of 
size L 3 where L = AAx with N equal to the number of cells in each direction and the grid spacing 



Ax = 0.2 fm. The fluid dynamic equations (18) are solved using the full (3+l)d SHarp And 



Smooth Transport Algorithm (SHASTA) ideal fluid dynamic code |644 165]. To ensure numerical 
stability we use a time step of At = 0.4 • Ax = 0.08 fm. 

Rewriting the equation of energy and momentum conservation for the coupled system with a 
source term S u for the fluid we obtain 

fyTf = -dp (Tr + Tf ) = S v . (23) 

After performing the fluid dynamical step for vanishing S u in the standard fashion, we subtract 
the sources from the energy and momentum density in the global rest frame of the fluid. This is 
especially interesting for an expanding medium and has been successfully implemented in previous 
studies of chiral fluid dynamics simulations without a Polyakov loop |37[ [38] . As will be shown 
later, this enables us to conserve the total energy for our Polyakov loop extended model, see section 

IIVBi 

To implement the source term in the numerical simulation we apply the same strategy as in 
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. Equations (19) and (20) provide us with the energy density dissipated into the heat bath. 

dQ q q „ \ „ f dQ, 



Ae di . 



da 



+ rj a d t cr d t a + 



gg 



Of 



+ m d t iT' d, 



At 



(24) 



The energy transfer due to stochastic fluctuations Aefl uc can be calculated by comparing this term 
to the numerically obtained energy difference in the sigma and Polyakov loop fields before and after 
each time step. The energy density of the sigma field is given by the sum of a potential, kinetic 
and fluctuation term while for the Polyakov loop it is associated only with the potential energy 
due to the lack of a kinetic term in the Lagrangian. 



TT/ , 1 [day 1 /-* >2 



ee = U{1) . 

Altogether this gives us the zero component of the source term 



1 

At 



(Ae d iss + Ae fl 

u 



(25) 
(26) 

(27) 



The spatial components accounting for momentum transfer are calculated in an analogous fashion. 
The momentum transfer due to dissipative processes is given by 



AM, 



diss 



+ r^cv) W + f^f + m d t tT 2 ) Yl 



At, 



(28) 



and AMfj uc is obtained by comparing this to the change in the momentum density of the fields 
which is vanishing for the Polyakov loop. 



M a = d t aVa 
M e = . 



This gives us the spatial part of the source term 

S 



^ (AMdi SS + AM fluc 



(29) 
(30) 



(31) 



The Langevin equations of motion for the fields a and i are solved using a staggered leap-frog 
algorithm, cf. [66]. For this algorithm the time step is chosen smaller than for the propagation of 
the fluid. We use a four times smaller At of 0.1 • Ax = 0.02 fm. I. e., four consecutive time steps 
of calculation for the fields are performed with the same fluid background, followed by one step for 
the fluid propagation. After that, the local temperatures T(x) are adjusted in each cell via a root 
finder of 

efluid(x) - e{a{x),£(x),T{x)} = . (32) 
These local temperatures then enter the local potentials ([7]), ( [To| ) and consequently the equations 



of motion ( 11 ) and ( 16 ) which are then used in the next time step for the propagation of the order 



parameter fields. 
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B. Results 



We are interested in the investigation of relaxational processes of the coupled system of fields 
and fluid after a temperature quench. Using a cubic 64 3 grid with periodic boundary conditions we 
set the temperature T[ n i uniformly to a value above T c and then initialize the sigma and Polyakov 
loop field at their equilibrium values including thermal fluctuations with a variance given by 



T 1 

V m 
T 1 

V m 



(^ 2 ) = T7^> (33) 



(34) 



Here we have defined the mass of the Polyakov loop rri£ analogously to the sigma mass as 

(35) 



•nip 



2 1 d 2 V cS 



e t 2 dt 2 

We then quench the temperature to various values T < T c and initialize the quark heat bath by 



calculating its energy density and pressure out of the given quantities T, a, I via equations (21) 



and (22). After that we let the system evolve according to equations (11), (16) and (23). Fields 
and fluid now influence each other in the following way: The amount of energy that the fields lose 
through damping gets transferred to the fluid which causes an adjustment of the temperature on 



account of equation ( 32 ) . This new temperature then reshapes the thermodynamic potential that 
influences the dynamics of the fields. In this kind of box calculations pressure gradients in the fluid 
are small, so we expect the dynamics to be governed by the fields. 

We are now interested in the relaxational behavior of the sigma field and the Polyakov loop 
comparing volume averages over all cells of the box as they evolve in time for different quench 
temperatures and two transition scenarios. The volume averages in a single event are defined as 

< CT > = a W ' W = £ »k , (36) 

ijk ijk 

where (lijk) is the instantaneous value of the sigma field (Polyakov loop) in a cell with coor- 
dinates i, j, k. These values are furthermore averaged over N e events with different noise configu- 
rations 

W = F &»- (Z) = W J2( £ )n- ( 37 ) 
n=l n=l 

Results for both transition scenarios are shown in fig. [3} They show the evolution of the noise 
and volume averaged value of the sigma field in time for different quenching temperatures T. 
The solid curves indicate the case where the system relaxes near the corresponding transition 
temperature. In both cases the dynamics is slowed down, however for different reasons. At the 
first order phase transition the barrier that separates minima near T c is responsible for a significant 
delay in the relaxation dynamics. For the second order transition we observe critical slowing down. 
This is inherent in the model due to the vanishing of the damping coefficient around the critical 
temperature that causes the system to oscillate around the equilibrium state and prolong the 
relaxation time up to infinity. 

The average values of the Polyakov loop as a function of time are shown in fig. [4} Here we find 
again a prolongation of the relaxation time near the first order phase transition where the average 
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FIG. 3. [3a] Equilibration of the sigma field for several quench temperature T < T c through the first order 
transition. The barrier between the minima in the potential increases the relaxation time when the system 
relaxes near T c = 172.9 MeV. We chose T ini = 180 MeV.|3b]Equilibrat ion of the sigma field for several quench 
temperature T < T c through the CP. Critical slowing down delays the dynamics and causes oscillations 
around the flat minimum when the system relaxes near T c — 180.5 MeV. We chose Tj ni = 186 MeV. 
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FIG. 4. 4a Equilibration of the Polyakov loop for several temperature quenches T < T c through the first 
order transition. We chose Ti n ; = 180 MeV. |4b] Equilibration of the Polyakov loop for several temperature 
quenches T <T C through the CP. We chose T ini = 186 MeV. 



value is slowly growing until t = 25 fm while for the other quenching temperatures, the system is 
equilibrated after 5 fm. For a CP scenario, we observe the same effect than for the sigma field: 
critical slowing down near the transition temperature, nevertheless with a small amplitude. 

We performed these simulations with various damping coefficients for the Polyakov loop r\i 
ranging from 1/fm up to 10/fm. A significant difference in the relaxational behavior could only be 
observed in the case where the system equilibrated near the first order phase transition where a 
larger value of rji caused a larger relaxation time. In all other cases the results are not sensitive to 
the choice of damping. Therefore we may consider our choice of r\i = 5/fm as justified, especially 
for the case of an expanding hot medium, which we finally aim to describe. 

The fluctuations of the order parameters can be analyzed by calculating their intensity N a and 
Ng. For the sigma field this quantity is given by [381 EH [68] 

dN„ a\a k = ul\5a k \ 2 + \d t a k \ 2 

d 3 k (27r) 3 2a; fc (27r) 3 2w fc ' { ' 

where at and a k are the Fourier coefficients of the expansion of the sigma field around its equilibrium 
value 5(7 = c — a eq and of the conjugate momentum field dto. The energy of the fc-th mode is 
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FIG. 5. 5a Intensity of sigma fluctuations during the transition process at t — 12 fm (first order) and t — 3 f m 
(CP). 5b Intensity of sigma fluctuations after equilibration at t = 24 fm. In the CP scenario we find an 



enhancement of the soft modes. [5c] Intensity of Polyakov loop fluctuations during the transition process at 
t = 12 fm (first order) and t = 3 fm (CP). 5d Intensity of Polyakov loop fluctuations after equilibration at 
t — 24 fm. In the CP scenario we find an enhancement of the soft modes. 



Uk = y k 2 + m 2 . We use an analogous definition for the Polyakov loop 



d 3 k 



^l\5£ k \ 2 + \d t i 
(27r) 3 2w fe 



(39) 



with ujk = \ k 2 + to 2 . , although this field formally has no kinetic energy term (see discussion at the 



end of section II C). In equilibrium the quantities N a , Nn may be interpreted as particle numbers, 
but to avoid confusion we call them 'intensities of fluctuations'. Several histograms of N a and 



Ng as a function of the wave number \k\ evaluated at different times are shown in figs. 5a, 5b for 



the sigma field and figs. 5c, 5d for the Polyakov loop. The figures on the left hand side show the 
intensity at early times during the transition process. We see that here fluctuations at the first 
order transition are clearly enhanced compared to the CP scenario. On the right hand side the 
intensities are shown for the time t = 24 fm after the system has equilibrated. Here we see that the 
long wavelength modes of both order parameters are enhanced at the CP, a typical and well-known 
critical phenomenon [16J. 



IV. DYNAMICS IN AN EXPANDING MEDIUM 



Here we are interested in the coupled dynamics of a system which is not confined in a box but 
freely expands into vacuum, similar to what happens after a heavy-ion collision. We study the 
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relaxational behavior of the sigma field and Polyakov loop during the nonequilibrium evolution 
and investigate the influence of energy-momentum exchange on the temperature evolution in both 
transition scenarios. 



A. Numerical implementation 

For this simulation we use a 128 3 grid and initialize in its center a droplet which is ellipsoidal in 
the x-y-plane and uniform in z-direction. This resembles the almond shape of the overlap region of 
two colliding nuclei. The droplet has a temperature of T[ Q { = 200 MeV, well above both transition 
temperatures and is smoothed by a Woods-Saxon distribution function at its edges: 

T- ■ 

T(x t = 0) = — (40) 

(1 + exp[(f - R)/a]) (1 + exp[(|z| - l z )/~a]) ' 



Here, f = y / x 2 + y 2 , a = 0.6 fm is the thickness of the transition layer to vacuum and 



"'" r ^ 



y/b?x 2 +a 2 y 2 ' _ (-42) 

a, f = 



R 



The ellipsoidal parameters are chosen as a = r& — b/2 and b = \J r 2 A — 6/4, = 6.5 fm denoting 

the radius of the two nuclei and b = 6 fm the impact parameter. The extent in z-direction is 
2l z = 12 fm. The sigma and Polyakov loop fields are initialized with their respective equilibrium 
distribution, then the energy density and pressure of the quark fluid are calculated. We choose an 
initial velocity profile with v z (x, t) = \z\/l z ■ v max with f max = 0.2. The transverse velocities v x and 
v y are set to zero in the beginning. 



B. Energy conservation 



We let the system expand by full 3 + 1 dimensional fluid dynamics and check the conservation 
of the total energy throughout the evolution. The energy of the sigma and Polyakov loop field are 
given by (25) and (26). Fig. [6] shows the total energy and the partial energies of the quark fluid, 
the sigma field and the Polyakov loop during the fluid dynamical expansion. For both scenarios 
the total energy is well conserved until the quark fluid reaches the edges of the computational grid 
after 8 fm. 



C. Supercooling and reheating 



During the fluid dynamic expansion we extract the average temperature (T), sigma field (a) 
and Polyakov loop (£) in a central cubic volume of lfm 3 inside the hot matter as a function of 



time. The results are shown in fig. 7a for the first order transition and in fig. 7b for a scenario 
with a CP. One can observe significant differences in the evolution of the average temperatures 



between both scenarios: For the case of the first order transition, fig. 7a, a reheating occurs after 
6 fm as a consequence of the formation of a supercooled phase below the transition temperature. 
We see that as the average temperature falls below T c , the average values of the sigma field and 
Polyakov loop remain close to their high temperature values around cr/ f w = 0.1 and I = 0.4. This 
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FIG. 6. [6a] The total energy as composed of the energy of the sigma field, Polyakov loop and the quark fluid 



of an expanding system for a first order transition scenario. 6b The total energy as composed of the energy 
of the sigma field, Polyakov loop and the quark fluid of an expanding system for a scenario with a CP. 



supercooled state decays after about 2 fm to the global minimum and transfers its energy into the 
fluid which consequently causes an increase in the average temperature. 

For the CP scenario, fig. |7b| no reheating effect is observed. The temperature decreases mono- 
tonically with only a small plateau well below the transition temperature where the dynamics 
slightly slows down due to the flat shape of the effective potential. The evolution of the averaged 
fields, especially for (a), proceeds less rapidly than at the first order phase transition. 
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FIG. 7. [7a|Evolution of the average temperature, sigma field and Polyakov loop in a system evolving through 
the first order transition. Supercooling followed by reheating can be observed. The horizontal line denotes 



the critical temperature. 7b Evolution of the average temperature, sigma field and Polyakov loop in a system 
evolving through the CP. The temperature decreases monotonically with a small plateau slightly below T c . 
The horizontal line denotes the critical temperature. 



D. Domain formation 



In our previous calculations we used the correlation functions (14) and (17) assuming that 
the stochastic noise fields in neighboring cells are not correlated as expressed by the spatial delta 
function. This leads to results that correspond to averaging over many events, see [38]. However, 
to better understand the differences between the two scenarios it is instructive to compare the 
microscopic structure of individual events. This requires a more consistent treatment of the field 
fluctuations by implementing realistic correlation lengths. Fig. [8] shows the correlation length of 
the sigma field £ = \jm a as a function of time for the expansion through both the CP and the 
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first order phase transition. Here, m a is calculated out of the volume averaged temperature in 
the previous section. The effect of such an averaging on the correlation length in inhomogeneous 
systems has been discussed in |69]. For the first order transition, £ lies in a range of 0.3 — 0.5 fm 
for the whole evolution while in the CP scenario it reaches a peak of about 1.5 fm when the system 
crosses the transition temperature after t = 3.2 fm, cf. fig. |7b| 

Below we present results obtained by correlating the noise fields over the spacelike distance 
£ = l/m a in each time step. The numerical procedure that implements such correlations has been 
described in [35J. For each spatial cell we perform an averaging of the randomly distributed noise 
field over a surrounding cube of linear size £ = nAx via 



C{x) = \ V £ a (x + iAxei + jAxe 2 + kAxe 3 ) . (42) 
i,j,k=0,...,n—l 



As this procedure weakens the fluctuations, i. e. (<5£^ 2 ) ^ we have to rescale the noise field 

in order to obtain again the initial distribution width: 



C(Z) = CW\l^l- (43) 



An analogous correlation procedure is done for the Polyakov loop, correlating & over the distance 
l/rri£. The following figures show spatial distributions in the transverse plane (z = 0), for the sigma 
field in fig. |9j the Polyakov loop in fig. 10 and for the energy density in fig. [TT] We chose a time 



of t = 4 fm for the first order scenario corresponding to the onset of the transition process where 



(T) = T c , cf. fig. 7a. Here we expect to observe phase coexistence, i. e. domains of the chirally 



broken phase in the chirally symmetric background or vice versa. For the CP scenario we chose 



a time of t = 3.2 fm, where again (T) = T c , cf. fig. 7b and furthermore the correlation length 
reaches its maximum value as shown in fig. [8} In this scenario there are no degenerate or metastable 
phases around T c but one single equilibrium state for each temperature. We therefore expect a 
more regular structure with no domains. 

By inspecting the figures we indeed find striking difference between the two transition scenarios. 
At the first order transition, both the fields and fluid evolve irregularly. This effect is best observed 
in the sigma field, fig. [9a] where one can see the expected domains of the chirally broken phase 
embedded in the chirally symmetric background, see for instance the region around (x,y) = (2, 2). 
Here, the system produces large fluctuations in small spatial regions that are able to overcome 
the potential barrier and create these characteristic structures. Nevertheless, the different phases 



17 



are connected smoothly, there are no sharp boundaries between the domains. This is a result of 
the Laplacian in the equation of motion which smoothens the gradients. Also in the Polyakov 
loop, fig. |10a and even more in the energy density, fig. 11a, we observe a bumpy and irregular 
structure. Regions of high energy density are embedded in a background of lower energy density 
in the periphery, like for instance around (x,y) = (—2,5). These embedded regions are then going 
to grow and merge, leaving small islands of the symmetric phase which gradually shrink until the 
whole system is in the low temperature chirally broken phase. 

On the other hand, the evolution through the CP proceeds smoothly, the regular ellipsoidal 
structure is preserved. Especially at the transition point, where the correlation length grows large, 



we see a smooth shape in both the fields and the energy density of the quark fluid, see figs. 9b, 10b 



and lib 



data, e. 



This striking difference in the event structure should manifest itself in the experimental 
;. in the non-statistical multiplicity fluctuations of produced hadrons [70] . 




V. SUMMARY AND OUTLOOK 

We presented a fully dynamical model to study the chiral and deconfinement phase transitions of 
QCD simultaneously. In our previous studies of chiral fluid dynamics we have derived the coupled 
dynamics of the sigma field and the quark heat bath within a self-consistent 2PI effective action 
approach. The results were now extended for the Polyakov loop so that both order parameters are 
propagated via Langevin equations taking into account the interaction with the quark heat bath 
via dissipation and noise. We assume that the structure of the source term for the Polyakov loop 
is analogous to that for the sigma field. During all simulations the total energy is well conserved. 

We studied the relaxational behavior of the coupled system for different quench scenarios in a 
box. For both the first order and the CP scenario, relaxation near the transition point is delayed. 
At the first order phase transition the relaxation process is significantly delayed due to the barrier 
in the thermodynamic potential. The transition actually starts only when this barrier disappears. 
Near the CP, when the sigma mass drops to zero and therefore the damping vanishes, we observe 
perpetual oscillations of the sigma field around the equilibrium value. These fluctuations are also 
visible in the Polyakov loop, although with a small amplitude. Performing a Fourier analysis of 
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the fluctuations for both order parameter fields we have observed another interesting peculiarity. 
While during the transition process the intensity of fluctuations is much stronger in a first order 
than in a CP scenario, the soft modes are stronger enhanced near the CP as compared with the 
first order transition if the systems are allowed to equilibrate. 

For the evolution of an expanding fluid through the first order transition we have found a 
clear evidence for the formation of a supercooled phase. Its decay later on leads to a substantial 
reheating of the quark fluid. That is in contrast to the simulation with a CP where the temper- 
ature decreases monotonically. Furthermore, during the onset of the first order transition, small 
domains of different phases coexist and create inhomogeneities in the energy density. In the CP 
scenario, where the correlation length grows large near the critical temperature, we observe a more 
homogeneous structure in the fields and fluid. A detailed analysis of domain formation at the first 
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order phase transition will be provided in future work. Furthermore, we will extend this model to 
finite baryochemical potential and study trajectories in the full T-/i-plane as well as fluctuations 
of baryon number densities. 

We think about several extensions and refinements of the model. The damping coefficient for 
the sigma field should include not only the a qq decay but also the decay into pions. The soft 
modes which are linked to the order parameter of the chiral transition are furthermore affected by 
interaction with the hard modes which would give another contribution to the heat bath and the 
dissipation process. For the Polyakov loop we used a purely phenomenological Langevin equation 
with a constant damping coefficient. This needs improvement, for instance by extracting this term 
from Monte Carlo simulations like it has been proposed in [32]. The present studies will allow for 
a better quantitative understanding of the signatures of the CP especially at FAIR energies. 
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